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We investigate the propagation of density- wave packets in a Bose-Hubbard model using the adap- 
tive time-dependent density-matrix renormalization group method. We discuss the decay of the 
amplitude with time and the dependence of the velocity on density, interaction strength and the 
height of the perturbation in a numerically exact way, covering arbitrary interactions and amplitudes 
of the perturbation. In addition, we investigate the effect of self-steepening due to the amplitude 
dependence of the velocity and discuss the possibilities for an experimental detection of the moving 
wave packet in time of flight pictures. By comparing the sound velocity to theoretical predictions, 
we determine the limits of a Gross-Pitaevskii or Bogoliubov type description and the regime where 
repulsive one-dimensional Bose gases exhibit fermionic behaviour. 



I. INTRODUCTION 

The study of strong interactions in one-dimensional 
Bose gases has recently attracted considerable interest, 
in particular the suggestion of Petrov et al. that 
in sufficiently dilute gases a regime appears in which 
ID bosons exhibit properties similar to that of a non- 
interacting Fermi gas. Following the realization of single 
mode atomic wires by using strong 2D optical lattices [3 , 
this so-called Tonks gas regime has indeed been seen in 
recent experiments U Q ■ 

Our aim in the present work is to study the prop- 
agation of density waves in strongly interacting one- 
dimensional Bose-Einstein condensates. Quite generally, 
the low-lying excitations in a Bose-Einstein condensate 
are sound-like and correspond to fluctuations of the con- 
densate phase j5j . The associated sound velocity depends 
on both the density and interaction strength and is diffi- 
cult to calculate microscopically in general. Beyond the 
weak interaction limit, where a Gross-Pitaevskii or Bo- 
goliubov description applies, very few results are avail- 
able, except for the particular case of one dimension. 
In that case an exact solution for the ground state and 
the elementary excitations is available for the continuum 
model with a short-range interaction through the well 
known Lieb-Liniger solution of the ID Bose gas 0, Q- 
Experimentally, density perturbations can be created by 
applying a localized potential to the system with a far de- 
tuned laser beam Alternatively, a phase imprinting 
method can be used, which allows to create solitonic ex- 
citations [TollTH. 

In our present work, we study the evolution of density- 
wave packets in a system of ultracold bosons which are 



subject to an optical lattice along the axial direction. In 
previous studies, the motion of Gaussian wave packets 
has been investigated theoretically for small density per- 
turbations or broad perturbations in three dimensions 
both with and without an optical lattice 0, 0, IT^ . 
These investigations were confined to the regime of weak 
interactions, describing properly systems with many par- 
ticles per site. Here, we focus on the case of one- 
dimensional systems at low filling, i.e. with approxi- 
mately one or less than one particle per site on average. 
This regime is of particular interest, since it allows one 
to study the behaviour of sound waves near the transi- 
tion from a superfluid to a Mott-insulating regime, as has 
been been realized experimentally by Stoferle et al. |l5j . 

As first pointed out by Jaksch et al., ultracold bosons 
in an optical lattice provide a perfect realization of 
the Bose-Hubbard model (Eq. which con- 

tains the interplay between their kinetic energy and their 
on-site repulsive interaction. The recently developed 
adaptive time-dependent density-matrix renormalization 
group method (adaptive t-DMRG) [H El [l^ is used 
to calculate the time-evolution of wave packets. This 
method allows us to find the time-evolution for both 
weak and strong coupling. In particular, it works best 
in an intermediate interaction regime, where other meth- 
ods are not reliable. We focus our investigation on the 
decay of the amplitude with time and on the sound ve- 
locity, i.e. the velocity of propagation of an infinitesimal 
perturbation. In addition, we determine the velocity of 
propagation of a perturbation with finite amplitude, thus 
entering nonlinear effects which are difficult to discuss 
analytically even in one dimension. 

We compare our numerical results in the limits of weak 
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and strong interaction to different approximations: For 
weak interactions a continuum description is applied, 
which leads to a system of bosons with (5-interaction, 
the Lieb-Liniger model Q . We compare the resulting 
sound velocity with our results and find good agreement 
up to intermediate interaction strength. A further sim- 
plification is obtained by treating the Lieb-Liniger model 
in a hydrodynamical approach. The sound velocity de- 
termined by this approach is that of a Gross-Pitaevskii 
type description. It agrees with our result only for rather 
small interaction strengths. In the limit of strong inter- 
actions and at low fillings, the Bose-Hubbard model can 
be mapped onto a model of spinless fermions ■ As 
expected, our numerical results for the sound velocity in 
this limit smoothly approach the value predicted from 
this mapping to fermions |2l| . 

The paper is organized as follows: we first introduce 
the Bose-Hubbard model, the analytical approximations 
and the numerical method used. Then we investigate 
the motion of the wave packet. We analyse the decay of 
the amplitude of the perturbation and the dependence of 
the velocity, in particular the sound velocity, on system 
parameters like the background density, the interaction 
strength and the height of the perturbation. Finally, we 
study how the presence of a wave packet can be detected 
experimentally from the interference pattern in a time of 
flight experiment. 



II. MODEL 

The Hamiltonian of the Bose-Hubbard (BH) model is 
given by 

L-i [/ ^ ^ 

= - J ^ + h.c. + ~ 1) + H 

(1) 

where L is the number of sites in the chain, bl and bj 
are the creation and annihilation operators on site j and 
nj = is the number operator [22|. In the Hmit of 
strong interactions, m ^ 1 with u := C//J, the atoms 
tend to localize. At integer filling p — N/L — 1,2..., 
where N is the total number of bosons, an incompress- 
ible Mott insulating phase with locked density arises once 
u is increased beyond a critical value [uc ~ 3.37 for p = 1 
according to p3 | in the thermodynamic limit). For weak 
interaction one finds a compressible superfluid phase. Ex- 
perimentally 0, 0, 01 , the parameter u can be varied 
over several orders of magnitude by changing the lattice 
depth. This allows one to tune through a superfluid- 
Mott-insulator transition, as first realized by Greiner et 
al. in a 3D optical lattice |23|. As mentioned before, it is 
possible to generate additional localized potentials using 
laser beams. These external potentials are modeled by 
the last term in equation In the following we use 
units in which the lattice spacing a — \, the hopping 
J = 1, and h = 1. This means that times are measured 



in units of h/ J and velocities in units of aJ/h. 

III. ANALYTICAL APPROXIMATIONS 

For weak interactions, or quite generally for a descrip- 
tion of the long wavelength properties of a noncommen- 
surate superfluid state, the continuum limit can be per- 
formed by taking Ja^ — const and a 0. In this limit 
the Bose-Hubbard model becomes equivalent to the Lieb- 
Liniger model 6, 7] 

^ /^"^ 21?'^^*^'^)'' + |(*t(a;))2(vI/(a:))^ (2) 

a bosonic model with (^-interaction of strength g. In this 
limit, the hopping parameter of the lattice model is re- 
lated to the mass M of the atoms by Ja^ — ^ and the 
interaction strength by U a — g. 

Starting from this continuum model and considering 
the interaction in the mean field approximation, the 
Gross-Pitaevskii equation can be derived . Within this 
approximation, the motion of density waves is described 
by the two coupled equations (written in dimensionless 
form) 

9p ^ djvp) ^ ^ 
dt dx 

Here p is the atomic density satisfying J dx p — 1 and 
V the velocity field. This is a good description for sys- 
tems in high dimensions or one-dimensional systems with 
many particles per site. Linearising the equations one re- 
covers the results of the hydrodynamical approach Q. 

We now turn to the opposite limit of strong interac- 
tions. For low densities p < 1 and strong interactions, the 
BH-Hamiltonian can be mapped to an effective model of 
spinless fermions with correlated hopping and attractive 
interactions 

-— E ("^+1 + "^-1) + O iiJ/U)') , (4) 

where {cj,c|,} = Sj,j', anticommuting otherwise, and 
fij = CjCj. Due to the correction O ((?7/J)~^), this map- 
ping is only valid for U / J ^ 1. 

IV. METHOD 

To study the evolution of a free wave packet in a ho- 
mogeneous system we apply the recently developed adap- 
tive t-DMRG 51 [13. The adaptive t-DMRG is a nu- 
merical method based on the well known static DMRG 
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[25I and the time-evolving block-decimation pro- 

cedure (TEBD) developed by Vidal The method 

describes the time-evolution of wave-functions in an es- 
sentially exact manner (for a detailed error analysis see 
[2^). In the calculation, the infinite-dimensional bosonic 
Hilbert space on a single site is truncated to a finite value 
Nb ■ We checked the consistency of our results by vary- 
ing Nb- For a chain of length L — 32 and not too high 
density, the results for Nb ~ 6 and Nb — 9 agreed well. 
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V. PREPARATION OF AN INITIAL STATE 

To prepare a density perturbation in our system we 
apply, for t < 0, an external potential Sj of Gaussian 
form, 

ej(t) = -2r)pe-^"/2-'0(-i) , (5) 

which is switched off for times t > 0. For weak perturba- 
tions, this potential creates an approximately Gaussian 
density packet 

Pj(i<0) = po(l + 2r;e-^"/(2'^')). (6) 

Note the difference between the parameters a and fj, 
which are used to describe the applied potential, and 
the parameters cr and 77, which determine the resulting 
density profile. For weak perturbations a = a, and 77 
is related to via the compressibility dp/dfi ^ 1/U. 
The background density po differs from the filling p not 
only by the effect of the perturbation but also by bound- 
ary effects. One constraint for the description of the 
time-evolution of a wave packet by the Bose-Hubbard 
model is that the bosons should not be excited to higher- 
lying energy bands induced by the periodic potential 
of the optical lattice. Hence it is valid as long as the 
additional energy by the perturbation is much smaller 
than the level spacing of the energy bands. The en- 
ergy change induced by the perturbation consists of two 
contributions: the change in the interaction energy and 
the change in the kinetic energy. The first can be ap- 
proximated by AE^int — 2pApU, with Ap ~ 7717 and 
U « i^Trash"^ /M) J (Px\w{x)\'^ , where as is the scatter- 
ing length and wlx) is the Wannier function. The ki- 
netic energy is dominated by the fast oscillations induced 
by the periodic lattice potential as long as the change 
in the density by the perturbation varies more slowly. 
Hence an upper bound for the change in the kinetic en- 
ergy is given by AEkin ^ JAp. In total we demand 
that AE ~ UAp{J/U + 2p) < hiy, where hiy is the en- 
ergy level spacing obtained approximating the wells by 
parabolic potentials. For p ^ 1 and J/U < 1, this con- 

dition is obeyed provided that rja <ti ^ 10, where 
a_L and a|| are the oscillator length perpendicular and 
parallel to the one-dimensional tubes. 



FIG. 1: Snapshots of the evolution of the density distribution 
are shown at different times. At f = 0, a Gaussian wave 
packet is present in the center of the system. It splits up 
into two packets which move with the same speed in opposite 
directions. 



VI. EVOLUTION OF THE WAVE PACKET 

A simple description of the evolution of a Gaussian 
wave packet for weak interactions can be obtained from 
a hydrodynamical approach. Linearizing equations (3), 
one obtains a linear wave equation. An initially Gaussian 
wave packet therefore shows a time evolution of the form: 

p(x,i) =po[l + ??(e-("-"*''/'"' -He-("+^*)'/2"')]. (7) 

The wave packet at < = thus splits into two packets, 
which travel with the same speed in opposite directions. 
Indeed this is the behaviour found in our simulations at 
weak coupling. Fig. shows snapshots of the evolution 
of a density- wave packet created at time t = 0. When 
the wave packets reach the boundaries, they are reflected 
back and after some time they meet again in the cen- 
ter of the system. The evolution of the density for up 
to four reflections is shown in Fig. |21by a density plot, 
i.e. the height of the density is encoded in a greyscale 
scheme. The bright lines indicate the motion of the wave 
packet, which splits into two packets moving towards the 
boundaries. After some time the pattern becomes less 
pronounced and a substructure arises due to the reflec- 
tion and scattering of the wave packets. 



VII. HEIGHT OF THE AMPLITUDE 

Damski has shown that, within the linearized 
Gross-Pitaevskii equation, the amplitude of the pertur- 
bation stays constant in time and equals po{l + i]). A 
decay of the amplitude in this approximation thus only 
occurs when nonlinearities become relevant. One of the 
origins of nonlinearity is the last term in Eq. Q , the so- 
called quantum pressure term. It arises from the kinetic 
energy term and describes a restoring force due to spa- 
tial variations in the magnitude of the wave function of 
the condensate. It becomes important if the length scale 
of spatial variations is of the order of the healing length 
^ = a/(-\/27po), where 7 is the dimensionless interaction 
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FIG. 2: Here the evolution of a density-wave packet is shown 
in a density plot. A linear grey scale is used, bright mean- 
ing higher densities. The bright lines correspond to the wave 
packets first splitting up moving towards the boundaries, be- 
ing refiected by the boundaries and meeting again in the cen- 
ter of the system, where the cycle starts again. After some 
reflections a substructure arises due to boundary effects and 
packet interactions. 



strength defined by 7 = . Hence a decay of narrow 
or high wave packets is expected even without an exter- 
nal potential. In agreement with this quaUtative picture, 
our numerical results for the Bose-Hubbard show that the 
decay becomes faster if (i) the width of the perturbation 
is narrower, and (ii) if the amplitude of the perturbation 
is higher. As an example, in Fig. O the decay of the 
amplitude is shown for different amplitude heights and 
widths. Both plots show a very rapid decrease for small 
times (in (a) for i < 1 and in (b) for t < 2), which is due 
to the splitting of the wave packet. For larger times after 
the two wave packets are separated, the decay is approx- 
imately linear in time (this might be just the first contri- 
bution of a more complicated decay). Clearly, the decay 
of the amplitude of the initially small height 77 « 0.1 and 
width fj « 1.4 [Fig. 01 (b)] is much slower than the decay 
of the amplitude of the initial height 77 « 0.3 and width 
CT w 1 [Fig. O (a)]. The oscillations seen in the curve 
stem from the discrete structure of the lattice, since we 
plot the maximum value of the lattice occupancies over 
all lattice sites (and not the maximum of an fitted con- 
tinuous curve which could lie between two lattice sites). 

Due to the rather slow decay of small amplitudes, we 
determine in the following the values to be used for po, 
77, and fj, by fitting the initial wave packet at < = 
to the form given by Eq. (|SJ). Such a fit is shown in 
Fig. n at i = 0. The error that results from assuming 
a time-independent amplitude rj is negligible for small 
amplitudes and broad widths of the perturbation. The 
uncertainties of the numerical results for the density (de- 
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FIG. 3: The typical decay of the amplitude of the perturba- 
tion. We plot (pj)max — Po, i-e. the difference between the 
largest discrete site occupancy and the background occupa- 
tion. The steep decrease for small times (up to t ~ 1 in (a) 
and t « 2 in (b)) corresponds to the splitting of the density- 
wave packet into two packets moving into opposite directions. 
The small oscillations in the curve stem from the discreteness 
of the underlying lattice. A linear fit is shown as a first ap- 
proximation. For the lower and broader amplitude (b) a much 
slower decay is seen as for the amplitude (a). 



termined by convergence checks in the number of DMRG 
states m, the allowed number of bosons Nb per site, and 
the Trotter time step At), and the errors made when 
reading off the parameters from the fit are much smaller 
than the size of the symbols used for data points in our 
plots (see for example Fig. 



VIII. SOUND VELOCITY 

To investigate the dependence of the sound velocity on 
the background density pa and the interaction strength 
u in the Bose-Hubbard model, we create two small den- 
sity perturbations with low amplitudes, a "bright" one, 
i.e. ?7 > 0, and a "grey" one, i.e. rj < 0, {\ri\ < 0.02) at 
approximately the same background densities. Since the 
sound velocity is the velocity for an infinitesimal density 
perturbation and we simulate the motion of perturba- 
tions with finite amplitude, we interpolate between the 
two results for the velocity of the perturbations ±r] lin- 
early (this will be justified later on, see section below). 
The velocity is determined from the propagation of the 
maximum or the minimum of the density perturbation for 
±77, respectively. In Fig. ^the sound velocity is plotted 
as a function of the interaction strength at fixed back- 
ground density po « 0.52 (The background density can 
not be fixed easily to a certain value, since it depends on 
the total number of particles, the boundary effects and 
the perturbation. In our calculations it deviates from po 
at most by 0.01.). 

Our numerical results will be compared with the the- 
oretical predictions from (i) a hydrodynamical approach 
or the linearized Gross-Pitaevskii equation, (ii) the Bo- 
goliubov approximation for the continuum gas by Lieb 
and Liniger, and (iii) the results of the mapping onto a 
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spinless fermion model. 

(i) The sound velocity determined by a hydrodynamical 
approach is given by 

v{p,9) = V9P- (8) 

Using the relations of the continuum limit, the corre- 
sponding velocity in the lattice is 

v{po,u) = 2pQ^/^, (9) 

with 7iat = u/2pa being the lattice analogy of the dimen- 
sionless interaction. 

(ii) As will be shown below, a much wider range of appli- 
cability than (i) is obtained from the results of Lieb and 
Liniger for the continuous bosonic model (Eq. ^ with 
(5-interaction. They found two distinct modes of excita- 
tions, the usual Bogoliubov mode and the Lieb mode, 
which is associated with solitary waves |2^. At low 
momenta the dispersion relations for both modes have 
the same slope, which means that they propagate at the 
same sound velocity. The expression for the sound ve- 
locity can be obtained from the thermodynamic relation 
mv"^ = pdpH, with /z as the chemical potential of the 
ground state, which is calculated within the Bogoliubov 
approximation. This results in 

(10) 

where vj = ^^^^ is the analog of the bare 'Fermi' velocity. 
In order to relate that to the Bose-Hubbard model, we use 
the expressions obtained from the continuum limit, i.e. 
7 — > 7iat and vj — > wjjat = 27rpo- Within the continuum 
model, the numerical calculation of the sound velocity 
by Lieb and Liniger shows that that expression (|10|l is 
quantitatively correct up to 7 ^ 10. By contrast the 
hydrodynamical result © is valid only up to 7 « 1. 

(iii) For strong interactions the sound velocity obtained 
by a mapping on a spinless fermion model is given by [2H 

- -fF (^1 - ;^(pocos7rpo)^ (11) 

where the Fermi velocity of the lattice model is = 
2 siuTrpg. 

In Fig. ^we compare these predictions to our numer- 
ical results. We see that for small interaction strength, 
M < 1, i.e. 7iat ^ 1 (note that for po — 0.52 u w 7iat), 
the curves obtained using (i) and (ii) agree well with our 
numerical results. Around 7iat ~ 1 the mean field pre- 
diction (i) starts to grow too fast, while the Bogoliubov 
approximation (ii) remains close to the numerical results 
up to intermediate interaction strength 7iat ~ 4. For even 
higher interaction strength also (ii) starts to differ signif- 
icantly from our numerial results. This means that the 
lattice model starts to deviate from the continuum model 
since (ii) was a very good approximation for the contin- 
uum model up to 7 « 10. A breakdown of the continuum 



limit in this regime is expected, since the healing length ^ 
becomes of the order of the lattice spacing a and thus the 
discreteness of the lattice becomes relevant. The sound 
velocity in the lattice model always remains lower than 
in the continuum model. For higher interaction strength 
the numerical results approach the asymptotic value of 
prediction (iii). Note, that the prediction (iii) is only 
expected to become valid for even stronger interactions 
than shown here, since it is an expansion in In Fig. 
El we see that our numerical results up to intermediate 
interaction strength show the dependence on the back- 
ground density predicted by Eq. IjlOl) . Deviations from 
the predicted form occur for 7iat > 2, depending on the 
particular set of parameters u and po- This dependence 
of the breakdown of the continuum limit (^ becomes of 
the order of a) is due to the fact that the healing length 
^ does not only depend on po E^nd u in the combination 
given by 7iat. Therefore the deviations at smaller values 
of u arise for larger background densities. Alternatively, 
this may be expressed in the form shown in Fig. |31 the 
breakdown of the continuum limit occurs for larger u at 
smaller 7iat. 

To summarize, we find that the sound velocity as a 
function of the interaction strength shows a crossover 
between Eq. (|10() . where Vg/po only dependends on the 
combination of po and u given by 7iat, to a saturation at a 
value given by Eq. (|ll|l . In fact, a completely analogous 
behaviour appears in the average kinetic energy of the 
particles, allowing to identify the Tonks regime for quasi 
ID tubes of bosons which are radially confined by a 2D 
optical lattice of increasing strength jj| . The breakdown 
of the prediction Eq. (|10|l is due to the discreteness of 
the lattice model and takes place if the healing length 
becomes of the order of the lattice spacing. 

The results presented above were obtained using chain 
lengths between L — 32 and L — 48 sites. Our numer- 
ical results for the time-evolution of the density profile 
are converged in the number of states kept in the re- 
duced space m (taken between to = 64 and to = 96), 
which means that the Trotter error dominates the total 
error j^^. The errors in observables are very small (of 
the order of 0.0001) for the Trotter time steps between 
At = 0.01 — 0.05 and can savely be neglected in compari- 
son to the uncertainties introduced by the determination 
of the sound velocity: For small interaction strength the 
velocity is relatively low and the movement over a long 
time can be fitted such that the accuracy of the results 
is of the order of ±0.01 before interpolation between zLrj. 
For higher interaction strength, the uncertainty in the re- 
sults for the velocity increases (approximately O(±0.05) 
for u — 6). This has two reasons: first, the velocity in- 
creases such that the end of the chain is reached in a 
rather short time. Moreover, oscillations in the density 
distribution induced by the finite size of the chain be- 
come more important and disturb the free evolution of 
the wave packets. 
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FIG. 4: The dependence of the sound velocity at constant 
background density po = 0.52 on the interaction strength is 
shown. Our numerical results (+) are compared to (i) the 
results Eq. @ of the hydrodynamical approach, (ii) the sound 
velocity determined by Lieb and Liniger Eq. IIIUI , and (iii) the 
results Eq. (Illll for strong interaction strength obtained by 
mapping onto spinless fermions. The results of Eq. Illll . i.e. 
(iii), should become applicable for even stronger interactions 
than the ones shown here. 
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FIG. 5: The dependence of the sound velocity on the inter- 
action strength and the background density is shown up to 
intermediate interaction strength. To confirm the prediction 
(ii) (Eq. 110^ we plot the ratio Vs/{2po) verus ^lat = u/(2po). 



FIG. 6: The dependence of the velocity on the height of the 
amplitude rj. The velocity is scaled by l{po) to remove its 
dependence on the background density. 

77, the dependence is approximately linear. It may be 
parametrized by 077 + 6 where a = 0.8 and b = 1.1. This 
linear dependence justifies the previously applied linear 
interpolation between ±77 for the determination of the 
sound velocity. 

As a consequence of the fact that the velocity increases 
monotonically with the amplitude of the perturbation, 
the wave can undergo self-steepening and shock wave 
formation can occur |l2l IT^ . One example where the 
phenomenon of self-steepening can be seen for a "bright" 
perturbation is shown in Fig. |3(a). It can be seen that 
the form of the density wave becomes very asymmetric. 
The front of the wave steepens and the back becomes 
more shallow. An additional dip arises at the front of 
the wave packet. This might stem from the discreteness 
of our system. In the case of a "grey" perturbation [Fig. 
El(b)], the asymmetry develops the other way round; the 
front becomes more shallow and at the same time the 
back of the wave steepens. It should be emphasized, how- 
ever, that the perturbations taken here are very narrow 
and high to have a strong effect. The BH model might 
not be quantitatively applicable to describe such pertur- 
bations. 



X. MOMENTUM DISTRIBUTION 



IX. SELF-STEEPENING 

In Fig. inithe dependence of the velocity on the height 
of the initial density-perturbation amplitude is shown. 
The data are taken at fixed interaction strength u = 1 
and different background densities po- The dependence 
of the velocity on the density is taken out by dividing 
by l{po) = \/2A)(1 - ^75^)^^^ ^^i'^g °™ knowledge 
from the previous results (cf. Eq. ^1 with 7 = 7iat = 
u/2pQ, and u = 1). We see that for small amplitudes 



Experimentally, one way of detecting the density per- 
turbation is to take time-of-flight images [Tsl . Theo- 
retically the interference pattern can be determined from 
the Fourier transform of the one-particle density matrix 

/(fc) = 1/iV ^ e'(^-^>'=(6]6,,) , 

neglecting its slowly varying envelope [s^- In a homo- 
geneous system without a density perturbation a sharp 
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FIG. 7: The evolution of a narrow density-wave packet is 
shown for various fixed times. The wave packets undergo self 
steepening and assumes an symmetric form. The lines are 
guides to the eye. 




FIG. 8: On the left the interference pattern is shown for two 
different times. At t = only one sharp interference peak 
at fc = exists. For times t > further peaks at finite 
momentum k and —k arise which correspond to the moving 
wave packets. Here only the region > is shown, exploiting 
a symmetry under k ~k. On the right the difference of the 
interference pattern for t = 5 and t = is shown. Here the 
errors are of the order of a few percent. 



XI. SUMMARY AND OUTLOOK 



interference peak appears at low interaction strength due 
to the long range order in the one-particle density ma- 
trix. If the interaction increases beyond the point where a 
Mott-insulating phase is present, this peak broadens and 
decreases. Finally, for very strong interaction only a dif- 
fuse pattern is left 'si'l . In the presence of a density- wave 
packet, we find that a second interference peak appears 
at a finite momentum. In Fig. |Hlwe show the difference 
between an interference pattern at t = 0, where the den- 
sity wave is still in the center, and a later point, where the 
wave packets travel through the system. The possibility 
to resolve the second peak in the experiments depends 
on the parameters of the system. Specifically, the peak 
shown in Fig. ISJa) was calculated for a high amplitude 
of the density perturbation. This ensures that the mean 
number of bosons contributing to the second peak in the 
interference pattern is a sufficiently large fraction of the 
total boson number. In Fig. |Sl(b) the difference between 
the pattern at t = 5 and i = is shown. 

In the experimental realizations a parabolic trapping 
potential is present in addition to the periodic lattice. 
As a result, the background density is no longer homo- 
geneous. Since the sound velocity depends on the back- 
ground density, we expect it to vary for weak interactions 
according to ((117)) and for strong interactions according 
to Only in the region where the trap varies slowly 

enough that the background density is almost constant, 
we expect the trap to have neglegible effect on the motion 
of the wave packet. 



To summarize, we investigated the motion of a wave 
packet in a Bose-Hubbard model which describes the dy- 
mamics of density perturbations in ultracold bosons in 
an optical lattice with a filling close to one particle per 
site. In the limit of weak interaction, 7^1, the mo- 
tion of relatively broad and small perturbations can be 
described by the hydrodynamical approach or the lin- 
earized Gross-Pitaevskii equation. For intermediate in- 
teraction strength, however, the mean-field description 
breaks down while the results obtained from the corre- 
sponding continuum Lieb-Liniger model remain valid in 
this regime (7 < 4). For strong interactions, we found 
that the sound velocity is well approximated by a map- 
ping onto a spinless fermionic model. In addition, we 
found a linear dependence of the velocity on the height 
of the amplitude. This gives rise to effects like self steep- 
ening and shock wave formation, in agreement with ana- 
lytical predictions. Finally, we have shown that a density 
wave may be detected experimentally as an additional 
peak in the interference pattern. 

Let us conclude by mentioning a few open questions: 
In the exact solution of the continuum model by Lieb 
and Liniger there are in fact two independent types of 
excitations: one of them exhibits a generalized Bogoli- 
ubov type dispersion, which is linear at small momenta 
and crosses over to a quadratic free particle behaviour 
at large momenta. The other one only exists in a fi- 
nite momentum range. It has been later identified as 
the solitar y w ave of the nonlinear Schrodinger equation 
in ID |29l l32|. As was shown by Lieb and Liniger, the 
velocity of the dark solitons for repulsive interactions is 
always smaller than the linear sound velocity, coinciding 
with the latter only in the limit of long wavelengths. Ex- 
perimentally, dark solitons have been observed in quasi 
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ID Bose-Einstein condensates, and have been identified 
by the fact that their velocity depends on the imposed 
phase gradient 0, 0| . In the case of a deep lattice po- 
tential, as is studied here, solitary waves are predicted 
to appear in the weak coupling regime u <^ 1 and for 
sufficiently wide density perturbations which can be de- 
scribed by the ID nonlinear Schrodinger equation. In 
addition, the presence of a lattice potential implies that 
atoms with momenta near a reciprocal lattice vector ac- 
quire a negative effective mass. This leads to the ex- 
istence of bright ga p solitons, a subject of considerable 
current interest [33 . Is^ Issj ] , in particular in connection 
with instabilities for strongly driven optical lattices [36j . 
In this paper we focused our investigations mainly on the 
case of perturbations with small momenta, for which the 
two modes cannot be distinguished by their velocity. It 
is an open question to which extent the density waves in 



our simulations, can be interpreted as solitary waves and 
in particular what happens to these stable excitations 
in the regime of strong coupling, where the nonlinear 
Schrodinger equation no longer applies. 
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